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SUMMARY 

An implicit finite -difference procedure is presented for solving the compressible 
three-dimensional boundary-layer equations. The method is second-order accurate, 
unconditionally stable (conditional stability for reverse cross flow), and efficient from the 
viewpoint of computer storage and processing time (eOOOOg storage and 0.002 second per 
■ nodal point on the CDC 6600 computer). The Reynolds stress terms are modeled by (1) a 
single -layer mixing length model and (2) a two-layer eddy viscosity model. These 
models, although simple in concept, accurately predicted the equilibrium turbulent flow 
for the conditions considered. Numerical results are compared with experimental wall 
and profile data for a cone at an angle of attack larger than the cone semiapex angle. 
These comparisons clearly indicate that the numerical procedure and turbulence models 
accurately predict the experimental data with as few as 21 nodal points in the plane nor- 
mal to the wall boundary. Research continues in the areas of convergence accelerator 
techniques (reduction of computer processing time), turbulence modeling, and extension 
of the computer code to general configurations (general geometry package development). 

INTRODUCTION 

A current design and analysis requirement of the aerospace industry is develop- 
ment of accurate and efficient numerical techniques and corresponding user-oriented 
computer codes for solving the compressible three-dimensional laminar, transitional, 
and turbulent boundary-layer equations for flows over general configurations. These 
codes would substantially reduce the cost and time currently required for the develop- 
ment of advanced aircraft through the substitution of numerical simulation for time- 
consuming and expensive experimental simulation and testing. The research focused on 
three-dimensional boundary-layer flows can be attributed to (1) experience gained over 
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the past decade in developing user-oriented codes for two-dimensional and axisymmetric 
flows, including the development of numerical procedures and simple but accurate mean- 
field turbulence models and (2) the increased availability of large -storage high-speed 
digital computer systems. 

Experience gained in two-dimensional and axisymmetric bovindary-layer flows 
indicates that computer codes for three-dimensional turbulent flows, if they are to be 
accepted as design/analysis tools by the aerospace industry, will require efficient and 
accurate numerical_methods with_suitable-turbulence-models for the Reynolds stresses. 
Numerical experimentation and detailed experimental turbulent boundary-layer research 
has resulted in the development and verification of mean -field turbulence models (eddy 
viscosity/mixing length) for two-dimensional boundary-layer flows which are. sufficiently 
accurate for application to a broad range of flow and boundary conditions. (See refs. 1 
to 4.) The numerical techniques developed for this class of flow can be directly applied, 
with minor modifications, to three-dimensional flows; however, zones of dependence and 
independence must be carefully treated for three-dimensional flows (see refs. 5 and 6). 
Although the extension of two-dimensional mean-field turbulence models to three- 
dimensional flows appears to be straightforward, numerical experimentation in which 
numerical results are compared with accurate three-dimensional profile and wall data is 
required before any evidence can be produced to prove or disprove this assumption - ^ ‘ 
(refs, 7 to 9), ^ 

If one assumes that accurate and efficient numerical procedures together with suf- 
ficiently realistic turbulence models can be developed on the basis of experience with 
two-dimensional flows, a number of problem areas still remain to be solved for general 
three-dimensional boundary-layer flows; these include: (1) selection and development of 
an optimum boundary-layer coordinate system; (2) development of general transforma- 
tions which will remove numerical problems associated with the generation of initial data 
planes, reduce the growth of the boundary layer in the computational region, and reduce 
the sensitivity of the numerical procedure to mesh-point distributions in the two spatial 
surface coordinates; and (3) availability and/or use of accurate three-dimensional invis- 
cid flow -fie Id solutions which are required for edge boundary condition specification. 

(See ref. 10.) Problems associated with (1) and (3) make it manditory that the boundary - 
layer codes be coupled with the three-dimensional inviscid flow -field codes to avoid 
excessive and time -consuming. data manipulation as well as to provide the possibility of 
accounting for displacement surface effects on numerical results. For general aerospace 
configurations the problem associated with obtaining accurate three-dimensional inviscid 
flow-field data may weH be the most difficult. However, substantial progress has been 
made in this particular area of research. Progress has also been made in the area of 
three-dimensional boundary-layer flows over the past few years (refs. 11 to 15),. A crit- 
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ical review of computational techniques for boundary -layer flows (two- and three- * 
dimensional) is presented in reference 16. 

In this paper a technique under development at the NASA Langley Research Center 
for solving the compressible three-dimensional laminar, transitional, and turbulent 
bovindary-layer equations is presented. The advantages and disadvantages of the impli- 
cit finite-difference procedure and Crocco transformation are discussed. The Reynolds 
stress terms are modeled by two mean-field (scalar invar i^t) models. Numerical ^ 
results are presented and compared with experimental data to determine the validity of 
the simple turbulence models and the accuracy of the nunierical procedure. 
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SYMBOLS 


damping term in turbulence model (eq. (12)) 

coefficients in boundary condition on shear equation (see eq. (25b)) 
coefficient of geometric progression for mesh-point distribution. 




skin-friction coefficient 

specific heat at constant pressure 

Van Driest damping factor (eq. (11)) 
incremental arc length 



= u/ue 


= v/ve 


= w/ue 


hj,h 2 ,h 3 metric coefficients (eq. (5)) 
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ki,k2,. . coefficients in turbulence models (eqs. (9) to (17)) 


•^ff 

L 


effective thermal conductivity (eq, (8)) 
reference length 

mixing lengths (eqs. (10) and (17)) 


Moe 

N 

Npr 


free -stream Mach number 

number of mesh points in the plane normal to the wall boundary 
Prandtl number 


Npr,t 

P 

^,L 


static turbulent Prandtl number 
pressure 

Reynolds number based on reference length 


T 


Tioo 



temperature 

free -stream velocity 

total velocity at edge of boundary layer 


u,v,w velocity component in the I-, rj-, and ? -direction, respectively 
xj,X 2 >X 3 physical coordinates 


Yi,Y 2,. . .,Y6 coefficients for derivative relations (eqs. (33) to (38)) 
a angle of attack 

®1>®2»* • '>“5 coefficients in standard equation (eq. (24)) 
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coefficients in equation (40) 

ratio of specific heats 

normal intermittency function (eq. (13)) 

streamwise intermittency function 

incremental mesh-point spacing in the 1 , 77 ,^ coordinates 
incompressible displacement thickness (eq. (19)) 
eddy viscosity 

transformed normal coordinate 
transformed cross -flow coordinate 
= T/Te 


Kj^,)C 2 geodesic curvatures (eq. ( 6 )) 




molecular viscosity 


^eff 


effective viscosity (eq. (7)) 



I 

I 

P 

T 

$ 

<t> 


transformed streamwise coordinate 

similarity parameter 

density 

total shear stress 
shear parameter (eq. (23)) 
circumferential angle 



w dummy variable 

o)s surface streamline angle 


Subscripts: 

e bovindary-layer edge quantity 

i,j,k grid point indices 


p,q,r dummy indices 

w wall value (C = 1) 


^,7} direction of quantity 


A bar over a sjmibol designates a dimensional quantity. 

GOVERNING EQUATIONS 


The governing equations are written in general form as follows (see fig. 1 for coor- 
dinate system; bar over a symbol designates a dimensional quantity): 

Continuity 


A(h2h3p5) + l^^ihgpv) + ^O^ihgPw) = 0 
^ -momentum 

hi H h2 h3 9C ^ ■ phi ph3 9?\h3 9?/ 

? 7 -momentum 

J.9V + X9V+J.9V + _ uvici = -4- ^ + 4- -9-/^ 9 h\ 

hi 9| hg 9?7 h3 9C ^ ^ phg ph3 aCyh3 Hj 


( 1 ) 


( 2 ) 


( 3 ) 


Energy 


.1. ££ + .v_ 9T ^ _ 
hi 9| hg 9^ h3 H 


1 

PCp 



+ J-iS\ 
h2 



/ 

, j . a ^eff af\ 
hs 9?yh3 9^ J 


( 4 ) 
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where hj, h 2 , and h 3 represent the metric coefficients for the incremental arc length 
ds; that is 

d §2 = (hi d^) 2 + (h 2 d7j)2 + (hg dC) 2 

The parameters ici and Kg are the geodesic 
and 77 = Constant, respectively; namely 

hihg 

_ } 

- 1 
^2 ” ‘hjhg 

The governing system is completed with the perfect gas equation of state and Sutherland's 
molecular viscosity law. 

Closure of equations (1) to (4) requires that the effective viscosity and ther- 

mal conductivity kgff be expressed ip terms of the dependent variables. These rela- 
tions are formulated as follows: 

iTe„ = M(i+|r) (7) 

and ^ 


(5) 

curvatures of the curves ^ = Constant 


( 6 ) 


^eff = 


I £ 

Npr\ 7T Npr^t 


( 8 ) 


where /I, T, Npj., and Npr^t represent the molecular viscosity, eddy viscosity, 
Prandtl number, and static turbulent Prandtl number, respectively. The streamwise 
intermittency function F (ref. 17) models the transitional region of flow and is a func- 
tion of 4 and 77 ; 0 ^ F ^ 1. In the present analysis the initiation and completion of the 
transitional flow process are empirically specified; however, correlation relations could 
be directly incorporated into the computer code. The eddy viscosity is assumed to be a 
scalar fimction independent of coordinate direction (refs. 8 and 18). The following sim- 
ple scalar invariant turbulence models are considered: 


Single -layer mixing length model: 
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where 


= ko tanhfe D 
*3,e \*^2 ^3,eJ 


D=l-«p(.S) 


y = 




■ T=(v*vr 

Two-layer eddy viscosity model: 
Inner law 


^Inner “ ^^2 

Ihg L 


+ gsj 


ill 1/2 


Outer law 


^outer “ ^4PUe,t®*y 

where 


<>2 = kjXgD 



- 2\1/2 
+ Ve j 

'^0 

1 (u2 + v2)^/^ 

Ue,t 


( 10 ) 

( 11 ) 

(13) 

(OgXgSXj^c) 

(*3,c<*3=5^3,e) 

(18) 

(19) 


(14) 

(15) 

(16) 
(17) 


The point where the inner and outer laws are matched X 3 ^ is obtained from the conti- 
nuity of eddy viscosity. For the results presented in this paper, kj, k 2 , ks, k 4 , ks, 
k 0 , and Npj. j were assigned values of 0.435, 0.09, 26, 0.0168, 5, 0.78, and 0.95, respec- 
tively. These represent the classical values generally accepted for equilibrium two- 
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dimensional boundary-layer flows (see refs. 3 and 4); however, note that although the' 
assigned values are sufficient over a broad range of flow and wall boundary conditions, 
modifications are required for certain classes of flow. (See refs. 4 and 10, for example.) 

TRANSFORMATION 

The use of physical coordinates introduces a number of problems for three- 
dimensional boundary-layer flows which can be circumvented by the intro^ction of a 
suitable transformation. If physical coordinates are used, three main problems are 
encountered as follows: (1) the numerical procedure and resultant solution are sensitive 
to the mesh -point distribution in the two surface spatial coordinates (Ax]^ and Ax 2); 

(2) the growth of the boundary layer in the streamwise (xj] and cross-flow (X2) directions 
requires the addition of nodal points in the X3 -direction as the solution prepresses (these 
two factors result in excessive computer processing time and/or computer code logic); 
and (3) initial data planes cannot be generated where the initial boundary-layer thickness 
is zero (for example, at the tip of a sharp body). Consequently, in the present procedure 
a transformation is introduced which avoids these problems and, in addition, minimizes 
the computer processing time and storage requirements. 

Equations (1) to (19) are first nondimensionalized (see ref, 11 for definition of non- 
dimensional variables), and a similarity -type transformation is introduced for the normal 
coordinate and velocity as follows: ’ 

h3 = i^-^l>3 . . (20) 

V^,L 

ooW (21) 

where Uoo is the reference velocity and for a sharp cone J = \/l. The metric coeffi- 
cients h2 and h3 are arbitrary functions of the coordinates. In order to cast the 
equations into Crocco-type form, the following function is defined: 

^ = (1_F)^/^ (22) 

where F = u/ug. The continuity and ^-momentum equations are combined to form the 
shear equation, where the shear parameter $ is defined by 

$ = --Lii/i + 1 r) (23) 

h3 T V M / 

Consequently; , F is replaced by $ as a new dependent variable, and H = w/u^ is 
imcoupled from the system. The governing system of equations reduces to three coupled 
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nonlinear partial differenUal equations in 0, G, and $ tc^ether with an explicit alge- 
braic relationship for H, The system assumes the following form: 


^ + anfii + 
9^2 ^ ^ 


"3 + “4 If + 


“sl^= 0 

^ dT] 


(24) 


where w represents 9, G, and and o!]^, q! 2 , 0 : 3 , a^, and a 5 are nonlinear 
coefficients. , 


.__The.bovmdary.conditions on equation (24) are as follows: r 


when C = 0 

0 = 1 G=1 4> = 0 


(25a) 


when C = 1 

0 = 0 w or' =fa,v)^ 

V9?/w 

G = 0 ) 



where aj and a 2 are functions of geometry and the inviscid edge conditions. 


(25b) 


The primary advantage of the Crocco-type transformation is that the solution 
domain is bounded between the definite limits 0 ^ 5 1. The only disadvantage of the 

transformation is that velocity overshoot in F is not allowed; that is, F, must increase 
monotonically from the specified wall value (slip at the surface can be specified) to unity 
at the edge boundary. Edge vorticity and streamline swallowing are not considered in 
this paper. • 


SOLUTION TECHNIQUE 


Equation (24) is solved in an Iterative mode with a marching implicit finite- 
difference technique suggested by Dwyer (ref. 19) and modified by Krause (ref. 20). The 
method is second-order accurate and unconditionally stable (conditional stability for 
reverse cross flow; see ref. 20). For turbulent flows a minimum of two mesh points in 
the C -plane must be located in the viscous sublayer; consequently, a variable mesh-point 
distribution is used. In the present study a geometric progression is assumed; that is; 
^^k+l/^^k “ ^ ~ 2> 3, . . ., N - 1 (ref. 3). Variable mesh-point distributions 

are also used in the and 7?-planes to minimize the computer processing time and stor- 

age requirements. A schematic of the difference molecule is presented in figure 2. 
Equation (24) is written at the point (i-l/2,j,k) and solved for the values of the dependent 
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variables 9, - G, and H at the point (i,j,k). Consequently, the partial derivatives of 
equation (24) at the point (i-l/2,j,k) are expressed as follows: 






3?^ 




1-1/2, j,k 


/3a> \ 

A-l/2,i,l 


_V9^A-i,j,k 


(26) 


(27) 


/^\ ^ ‘^i,i,k - ‘^i-l,j,k 

\9^A-l/2,j,k ^^i-1 


(28) 


‘^i-l,j+l,k - ‘^i -l,j,k ^ ^i,j,k - ^i,j-l,k 


/3a)\ 

v^'^A-i/2,j,k 


Am 


Atj 


1-1 


(29) 


WS _1 4 Jc + OJs ,• If 

The derivative quantities in equations (26) and (27) are obtained from 

= YiUJp^q^r+1 - Y2Wp^q^r + Y30)p,q,r-1 

^P,q,r 


(30) 


(31) 


(If) " '^4‘^,q,r+l - YsWp^q,r - Ye‘^p,q,r-1 

\ ^/p,q,r 


(32) 


where 


Yl = 


Y2 = 


Y3 = 


ACr(ACr + A?r_l) 


^^r ^?r-l 


^?r-l(^^r + ^?r-l) 


(33) 

(34) 

(35) 


Y4 = 


^?r-l 




(36) 


29 



Y _ ACr-i -A^^ „ 




A^r_l(A?r + ^?r.l) 



For equally spaced mesh-point distributions in the ^ -plane, equations (33) to (38) assume 
the values 

Y3 = Yi 

Y5-.O „ • Y6 = Y4 

When a converged solution cannot be obtained at the most leeward plane, 0 = 180° (for 
example, for separation on leeward surface), a cubic Crank-Nicolson differencing scheme 
is used at the maximum 77 -station (ref. 14). If this procedure were not incorporated into 
the program logic, one 77 -station would be lost for each incremental because of the 
Krause differencing scheme; that is, equation (29) assumes the existence of a converged 
solution at the point (i-l,j+l,k) for K = 2, 3, , . N - 1, 

^ , The marching procedure cannot be initiated without the existence of two orthogonal 
initial data planes. For a sharp right circular cone these planes of initial data are gen- 
erated directly from the governing equations by using a second-order Crank-Nicolson 
scheme for the two planes ^ = 0, 0 = 77 i Vmax 0 i | ^ ^max’ V = ^ where simi- 

larity exists. A discussion of problems associated with obtaining initial data planes for 
general configurations is presented in reference 10. 

Substitution of equations (26) to (30) into equation (24) results in a system of coupled 
algebraic equations whose coefficient matrix is of tridiagonal form which can be effi- 
ciently solved for the dependent variables (Thomas' algorithm). The primary problem 
associated with equation (24) is that the coefficients (o!i, «2> etc.) are hignly nonlinear. 
The shear equation controls the convergence rate of the numerical procedure (iterations 
required) as the system is’ sequentially iterated. Equation (24) can be written for $ as 



9^2 $ 2 ja ?-^2 $ ^2 9 ? ^2 977 


where Pi, - P2t functions of geometry, inviscid edge conditions, and 

previous iterate values of the dependent variables F and 0 and their derivatives. The 
problem is further complicated by the inclusion of the turbulence models (eqs. (9) to (19)), 
since in the transformed plane ^ appears explicitly in the transformed relationships. 
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Consequently, the coefficients 03j, ^2f also depend on ^ for turbulent flows (for 
laminar flows this dependence is removed). The system of equations will not converge if 
the shear equation is written as shown in equation (40) because of the term. Con- 
vergence can be achieved by using a Taylor's series expansion of about the previ- 

ous iterate value $q; that is - 

i = ~y+o($2) (41) 

$ $G\ 9gI 

Substitution of equation (41) into equation (40) yields 



^2 , % 3^ , <^4 3^ - Q 
$2 3 ? ^2 an 


(42) 


Equation (42) converges in an average of five to seven iterations for high Reynolds num- 
ber turbulent flow. The wall boundary condition on $ (see eq. (25b)) also presents a 
problem since is unknown; however, the wall derivative relationship can be directly 
incorporated in the iterative solution procedure. In principle, it should be possible to 
reduce the average number of iterations substantially to a maximum of three, Itesearch 
continues in the areas of (1) restructuring equation (42), (2) treatment of the $ wall 
boundary conditions, and (3) the problem associated with ^ in the transformed turbu- 
lence models. Note however that the present procedure requires essentially the same 
processing time per mesh point (0,002 sec) as the Cebeci -Keller Box method (ref. 10) and 
that this time may be substantially reduced through convergence accelerator procedures 
and/or the inclusion of Newton -Rap hson iteration. 


RESULTS AND DISCUSSION 

The munerical procedure and turbulence models have been applied to a number of 
flows (current geometry limited to sharp right circular and elliptic cones). In this paper, 
numerical results are compared with experimental wall and profile data for a cone with a 
12.5° semiapex angle at an angle of attack of 15.75°, The free-stream Mach number, 
total pressure, and total temperature were 1.8, 172,4 kN/m2, and 294 K, respectively. 
Transition was assumed to be initiated and completed in the region 0.03 ^ ? 0.08 

(L = 105,6 cm). The adiabatic wall boundary condition was imposed on the energy equa- 
tion (see eq. (25b)); that is, (^] =0. No experimental data were input into the viscous 

flow solution. The inviscid pressure distribution Pg = Pg(|,7/) was obtained from a 
numerical solution of the three-dimensional inviscid flow equations. Experimental data 
for verification of the accuracy of the numerical procedure and turbulence models were 
obtained from reference 21. 
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The numerical results for F, G, and 9 are compared with experimental data in 
figure 3 for circumferential locations of <f> = 0°, 45°, 90°, and 135°. In order to eval- 
uate the effect of nodal -point spacing in the ?-plane, a parametric study was made for 
N = 301, 201, 101, 61, and 21 with ^?k+l/^?k ” results for N = 301 and 201 

were essentially identical, and those for N = 101 were within 0. 5 percent of the N - 301 
results. The agreement between the experimental and numerical results is very good for 
301 points and, in general, good for 21 points. The two turbulence models (eqs. (9) to (19)) 
produced essentially identical results, ^he two^layerjnodel^^esults preseMed in Jigure A- 
are for N = 301; however, the two-layer results for N = 21 were essentially identical 
to the N = 21 results of the single -layer model. A comparison of the numerical results 
for Cf^e presented in figure 4(a) indicates that the difference between the results for 
61 points and 301 points is approximately 1 percent and between 21 points and 301 points 
is approximately 3 percent. In figure 4(b), Cf^e is presented as a function of Fig- 
ures 3 and 4 indicate that as few as 21 points normal to the wall boundary can be used to 
obtain results to within 3 -percent accuracy (compared with N = 301 results). Numeri- 
cal results for surface streamline direction Wg = tan-^^Tjji^T^j obtained for N = 301 

and 21 are compared with experimental data in figure 5. The agreement is good consi- 
dering that the inviscid pressure distribution was obtained from the inviscid equations 
and not from experimental data; that is, displacement surface effects are not included in 
the viscous/inviscid calculations. 

The major points which should be noted in these comparisons are (1) that the numer- 
ical procedure is efficient and accurate and (2) that the turbulence models are satisfac- 
tory for high Reynolds number equilibrium turbulent boundary-layer flows. The Crocco- 
type transformation and the numerical procedure allow the generation of accurate solutions 
for a minimum of 21 points normal to the wall boundary. The computer code requires 
600008 storage (the i-l,j,k data plane is stored on disk) and approximately 0.002 sec- 
ond per grid point processing time on a CDC 6600 computer system. Current studies 
indicate that it may be possible to substantially reduce the processing time through con- 
vergence accelerators for the shear equation (eq. (42)) and/or the inclusion of a Newton- 
Raphson iteration procedure. The current program is comparable in both storage and 
processing time with the Cebeci-Keller Box method (ref. 10), 

CONCLUDING REMARKS 

Solutions of the compressible three-dimensional turbulent boundary-layer equations 
have been obtained and compared with experimental data. The agreement between the 
numerical results and experimental data indicates that accurate results can be obtained 
with a minimum of 21 nodal points in the plane normal to the wall boundary layer for high 
Reynolds number equilibrium turbulent flows. The turbulence models, although simple in 
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concept, were adequate for the class of flow considered; however, previous experience 
indicates that caution should be exercised in extending these models to more demanding 
boundary-layer flows. The numerical procedure is second-order'accurate and uncondi- 
tionally stable (conditional stability for reverse cross flow). The computer code requires 
60000g storage and approximately 0.002 second per ^id point processing time (CDC 6600 
computer system). Studies indicate that the processing time may be further reduced 
through convergence accelerator and Newton -Raphson iteration procedures; however, the 
current computer requirements of storage and speed compare favorably with other pro- 
cedures imder development. 
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Numerical results N 
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(a) Details of figure 3. 



Figure 3.- Profile comparisons. Moo = 1-8; Roo l = 25 x 10®; y = 1.4; 
z » 2.54 cm; L = 105.6 cm; Qg “ ^Ug^ + Vg^. 
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Figure 5.- Surface flow direction, xj/l = 1. 


